Simulations fitness data
Laure Olazcuaga, Nicolas Rode
25-February-2021
1 COLORS plot
colors <- c("LMER" = "#DCA237", "Anova" = "#459B76")
lines <- c("Balanced" = "solid", "Unbalanced" = "dotted")2 ESTIMATES VARIANCE REAL DATA
#############################
#PERFORMANCE PARAMETERS VALUE
#############################
data_PERF <- read.table(file=here::here("data", "DATACOMPLET_PERF.csv"), sep=",", header=TRUE)
data_PERF$Original_environment<-as.factor(data_PERF$Original_environment)
data_PERF <- data_PERF[data_PERF$Original_environment!="WT3",]
data_PERF <- data_PERF[data_PERF$Test_environment!="GF",]
data_PERF$Test_environment <- as.factor(data_PERF$Test_environment)
data_PERF$Population <- as.factor(data_PERF$Population)
data_PERF <- data_PERF[data_PERF$Test_environment!="Grape",]
data_PERF <- droplevels(data_PERF)
data_PERF$IndicG0<-as.numeric(ifelse (as.character(data_PERF$Gen) == "G0", 1, 0))
nfruit_data <- nlevels(data_PERF$Original_environment)
nhab_data <- nlevels(data_PERF$Test_environment)
npop_per_fruit_data <- nlevels(data_PERF$Population)/nfruit_data
nrep_data <- mean(tapply(data_PERF$Nb_adults,list(data_PERF$Population,
data_PERF$Test_environment,
data_PERF$Generation),length),na.rm = TRUE)
ntrial_data <- mean(data_PERF$Nb_eggs, na.rm = TRUE)
#POISSON VARIANCE VALUE
mperf_poisson <- lme4::lmer(log(Nb_adults+1) ~ 1+ (1|Population) +
(1|Original_environment:Test_environment) +
(1|Original_environment:Test_environment:IndicG0),
data = data_PERF)
sdpop_datapoisson <- sqrt(unlist(lme4::VarCorr(mperf_poisson))[1])
sdfruithab_datapoisson <- sqrt(unlist(lme4::VarCorr(mperf_poisson))[3])
sdfruithab_ng_datapoisson <- sqrt(unlist(lme4::VarCorr(mperf_poisson))[2])
sdfruithab_datapoisson <- sdfruithab_ng_datapoisson
sigma_datapoisson <- sigma(mperf_poisson)
#BINOMIAL VARIANCE VALUE
data_PERF_rate <- data_PERF[data_PERF$Nb_eggs>=data_PERF$Nb_adults,]
mperf_binomial <- lme4::lmer(asin(sqrt(Nb_adults/Nb_eggs)) ~ 1+ (1|Population) +
(1|Original_environment:Test_environment) +
(1|Original_environment:Test_environment:IndicG0),
data = data_PERF_rate)
sdpop_databinomial <- sqrt(unlist(lme4::VarCorr(mperf_binomial))[1])
sdfruithab_databinomial <- sqrt(unlist(lme4::VarCorr(mperf_binomial))[3])
sdfruithab_ng_databinomial <- sqrt(unlist(lme4::VarCorr(mperf_binomial))[2])
sdfruithab_databinomial <- sdfruithab_ng_databinomial
sigma_databinomial <- sigma(mperf_binomial)
###########################
#PREFERENCE PARAMETERS VALUE
###########################
data_PREF <- read.table(file=here::here("data", "DATACOMPLET_PREF.csv"), sep=",", header=TRUE)
data_PREF$Original_environment<-as.factor(data_PREF$Original_environment)
data_PREF <- data_PREF[data_PREF$Original_environment!="WT3",]
dim(data_PREF)## [1] 4428 12data_PREF <- data_PREF[data_PREF$Test_environment=="Cranberry"|
data_PREF$Test_environment=="Cherry"|
data_PREF$Test_environment=="Strawberry",]
data_PREF$Test_environment <- as.factor(data_PREF$Test_environment)
data_PREF$Population <- as.factor(data_PREF$Population)
data_PREF$BoxID <- as.factor(data_PREF$BoxID)
data_PREF$Nb_eggs <- as.numeric(as.character(data_PREF$Nb_eggs))
data_PREF <- droplevels(data_PREF)
data_PREF$IndicG0<-as.numeric(ifelse (as.character(data_PREF$Gen) == "G0", 1, 0))
nfruit_data_box <- nlevels(data_PREF$Original_environment)
nhab_data_box <- nlevels(data_PREF$Test_environment)
npop_per_fruit_data_box <- nlevels(data_PREF$Population)/nfruit_data
nrep_data_box <- mean(tapply(data_PREF$Nb_eggs,list(data_PREF$Population,
data_PREF$Test_environment,
data_PREF$Generation),length),na.rm = TRUE)
#POISSON BOX VARIANCE VALUE
mperf_poisson_box <- lme4::lmer(log(Nb_eggs+1) ~ 1+ (1|Population) +
(1|Original_environment:Test_environment) +
(1|Original_environment:Test_environment:IndicG0) +
(1|BoxID),
data = data_PREF)
sdbox_datapoisson_box <- sqrt(unlist(lme4::VarCorr(mperf_poisson_box))[1])
sdpop_datapoisson_box <- sqrt(unlist(lme4::VarCorr(mperf_poisson_box))[2])
sdfruithab_datapoisson_box <- sqrt(unlist(lme4::VarCorr(mperf_poisson_box))[4])
sdfruithab_ng_datapoisson_box <- sqrt(unlist(lme4::VarCorr(mperf_poisson_box))[3])
sdfruithab_datapoisson_box <- sdfruithab_ng_datapoisson_box
sigma_datapoisson_box <- sigma(mperf_poisson_box)#NUMBER OF SIMULS
# Power Analyses
nb_simul <- 5000
#Save the same number of simuls for each SA value
nb_per_SA <- 500
#
# #Number of simuls for each SA non-genetic value
# tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)
#False positive rate
nb_simul_fpr <- 1003 BALANCED AND UNBALANCED DESIGN
3.1 Normal distribution
3.1.1 Balanced
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", design = "balanced",
rho = -0.05, rho_ng = -0.05,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim1 <- as.data.frame(t(sim1))
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", design = "balanced",
rho = -0.5, rho_ng = -0.5,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim2 <- as.data.frame(t(sim2))
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal", design = "balanced",
rho = 0, rho_ng = 0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
#Add indic SA
sim <- add_indic_sign(sim)
#Add number of simul for each SA value to keep only 100 simuls per SA value
sim <- add_sim_number_SA(sim)
## Local adaptation (genetic)
tapply(sim$IndicGen[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], length)## -1.4 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2
## 1 3 6 13 16 44 69 100 100 100 100 100 100 100 100 100
## 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
## 100 100 100 100 100 100 100 100 100 100 100 61 26 13 5tapply(sim$IndicGen_aov[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], mean)## -1.4 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0300000 0.3900000
## 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## 0.5700000 0.5300000 0.6800000 0.6500000 0.6500000 0.6900000 0.6500000 0.7900000
## 1.1 1.2 1.3 1.4 1.5 1.6 1.7
## 0.8000000 0.8600000 0.8300000 0.8688525 0.8846154 0.9230769 0.8000000## Local adaptation (non-genetic)
tapply(sim$IndicNonGen[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)## -1.5 -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
## 0.1400000 0.1600000 0.2300000 0.2700000 0.3700000 0.3200000 0.5500000 0.5200000
## 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
## 0.5100000 0.6700000 0.6600000 0.6500000 0.7714286 0.5806452 0.8125000 1.0000000
## 1.8
## 0.6666667tapply(sim$IndicNonGen_aov[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)## -1.5 -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
## 0.1800000 0.1700000 0.2300000 0.2800000 0.3900000 0.3100000 0.5300000 0.5000000
## 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
## 0.5000000 0.6700000 0.6500000 0.6400000 0.7714286 0.5806452 0.8125000 1.0000000
## 1.8
## 0.6666667tapply(sim$IndicGen, sim$SA_Gen, length)## -1.4 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2
## 1 3 6 13 16 44 69 134 187 315 449 495 520 521 537 560
## 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
## 750 1163 1538 1756 1723 1447 1058 732 477 242 139 61 26 13 5tapply(sim$IndicNonGen, sim$SA_NonGen, length)## -1.5 -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1
## 1 1 1 3 15 18 34 75 147 218 292 405 481 530 570 527
## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
## 533 732 1126 1611 1840 1724 1400 1080 658 486 215 152 70 31 16 5
## 1.8
## 3####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)## -1.4 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2
## 1 3 6 13 16 44 69 134 187 315 449 495 520 521 537 560
## 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7
## 750 1163 1538 1756 1723 1447 1058 732 477 242 139 61 26 13 5sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
data_prop$design <- 0
data_prop_gen <- data_prop
#Plot
plot_genetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Normal data") +
theme_LO_sober +
ggtitle("Normal data") +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_genetic_normal
#################################
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 0
data_prop_nongen <- data_prop
#Plot
plot_nongenetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
theme_LO_sober +
ggtitle("Normal data") +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_nongenetic_normal
3.1.2 Unbalanced (low)
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal",
design = "unbalanced", disp_design = 1,
rho = -0.5, rho_ng = -0.5,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal",
design = "unbalanced", disp_design = 1,
rho = 0, rho_ng = 0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal",
design = "unbalanced", disp_design = 1,
rho = -0.05, rho_ng = -0.05,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA
sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
data_prop$design <- 1
data_prop_gen <- rbind(data_prop_gen,data_prop)
#Plot
plot_genetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Normal data") +
theme_LO_sober +
ggtitle("Normal data") +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_genetic_normal
#################################
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 1
data_prop_nongen <- rbind(data_prop_nongen,data_prop)
#Plot
plot_nongenetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
theme_LO_sober +
ggtitle("Normal data") +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_nongenetic_normal
3.1.3 Unbalanced (high)
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal",
design = "unbalanced", disp_design = 4,
rho = -0.5, rho_ng = -0.5,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal",
design = "unbalanced", disp_design = 4,
rho = 0, rho_ng = 0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "normal",
design = "unbalanced", disp_design = 4,
rho = -0.05, rho_ng = -0.05,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA
sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
data_prop$design <- 4
data_prop_gen <- rbind(data_prop_gen,data_prop)
#Plot
plot_genetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Normal data") +
theme_LO_sober +
ggtitle("Normal data") +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_genetic_normal
#################################
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 4
data_prop_nongen <- rbind(data_prop_nongen,data_prop)
#Plot
plot_nongenetic_normal<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
theme_LO_sober +
ggtitle("Normal data") +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_nongenetic_normal
3.1.4 Plot
data_proportion_gen<-tidyr::gather(data_prop_gen, "method", "proportion", 2:3)
data_proportion_gen$method <- plyr::revalue(data_proportion_gen$method,
c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_gen$design <- as.factor(data_proportion_gen$design)
#Plot genetic
plot_genetic_normal_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_gen,
aes(x = val_SA_genet, y = proportion, color = method, linetype = design)) +
geom_point(size = 2) +
geom_line(size = 1) +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Normal data") +
theme_LO_sober +
ggtitle("Normal data") +
labs(color = "Methods", linetype = "Unbalanced\ndesign") +
scale_color_manual(values = c("#459B76","#DCA237")) +
scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_genetic_normal_balancedvsunbalanced
#Data non-genetic
data_proportion_nongen<-tidyr::gather(data_prop_nongen, "method", "proportion", 2:3)
data_proportion_nongen$method <- plyr::revalue(data_proportion_nongen$method,
c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_nongen$design <- as.factor(data_proportion_nongen$design)
#Plot non-genetic
plot_nongenetic_normal_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_nongen,
aes(x = val_SA_nongenet, y = proportion, color = method, linetype = design)) +
geom_point(size = 2) +
geom_line(size = 1) +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
ggtitle("Normal data") +
theme_LO_sober +
ggtitle("Normal data") +
labs(color = "Methods", linetype = "Unbalanced\ndesign") +
scale_color_manual(values = c("#459B76","#DCA237")) +
scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_nongenetic_normal_balancedvsunbalanced
3.1.5 False positive rate analyses
3.1.5.1 No genetic and no non-genetic effects
sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "normal", design = "balanced",
rho=0, rho_ng=0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
## Local adaptation (genetic)
mean(sim$IndicGen)## [1] 0.02mean(sim$IndicGen_aov)## [1] 0.02## Local adaptation (non-genetic)
mean(sim$IndicNonGen)## [1] 0.01mean(sim$IndicNonGen_aov)## [1] 0.023.1.5.2 Genetic effects and no non-genetic effects
sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "normal", design = "balanced",
rho=-0.1, rho_ng=0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
## Local adaptation (genetic)
tapply(sim$IndicGen, sim$SA_Gen, length)## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2
## 3 4 14 6 18 17 13 12 7 4 2tapply(sim$IndicGen_aov, sim$SA_Gen, mean)## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
## 0.0000000 0.2500000 0.0000000 0.5000000 0.2222222 0.3529412 0.4615385 0.3333333
## 1 1.1 1.2
## 0.5714286 0.7500000 1.0000000## Local adaptation (non-genetic)
mean(sim$IndicNonGen)## [1] 0.02mean(sim$IndicNonGen_aov)## [1] 0.02tapply(sim$IndicNonGen, sim$SA_Gen, mean)## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
## 0.0000000 0.0000000 0.0000000 0.1666667 0.0000000 0.0000000 0.0000000 0.0000000
## 1 1.1 1.2
## 0.0000000 0.2500000 0.0000000tapply(sim$IndicNonGen_aov, sim$SA_Gen, mean)## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
## 0.0000000 0.0000000 0.0000000 0.1666667 0.0000000 0.0000000 0.0000000 0.0000000
## 1 1.1 1.2
## 0.1428571 0.0000000 0.0000000tapply(sim$IndicNonGen, sim$SA_NonGen, mean)## -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 0.1 0.2 0.3 0.4 0.5 0.6 0.7
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.6666667 0.0000000tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)## -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 0.1 0.2 0.3 0.4 0.5 0.6 0.7
## 0.0000000 0.1111111 0.0000000 0.0000000 0.0000000 0.3333333 0.00000003.1.5.3 Non-genetic effects and no genetic effects
sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "normal",design = "balanced",
rho=0, rho_ng=-0.1,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
## Local adaptation (genetic)
mean(sim$IndicGen)## [1] 0.1mean(sim$IndicGen_aov)## [1] 0.1tapply(sim$IndicGen, sim$SA_Gen, mean)## -1 -0.9 -0.7 -0.6 -0.5 -0.4 -0.3
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## -0.2 -0.1 0 0.1 0.2 0.3 0.4
## 0.00000000 0.00000000 0.00000000 0.08333333 0.00000000 0.12500000 0.20000000
## 0.5 0.6 0.7 0.9 1.2
## 0.33333333 0.50000000 0.33333333 0.50000000 1.00000000tapply(sim$IndicGen_aov, sim$SA_Gen, mean)## -1 -0.9 -0.7 -0.6 -0.5 -0.4 -0.3
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## -0.2 -0.1 0 0.1 0.2 0.3 0.4
## 0.00000000 0.00000000 0.00000000 0.08333333 0.00000000 0.12500000 0.20000000
## 0.5 0.6 0.7 0.9 1.2
## 0.33333333 0.50000000 0.33333333 0.50000000 1.00000000tapply(sim$IndicGen, sim$SA_NonGen, mean)## 0.2 0.3 0.4 0.5 0.6 0.7 0.8
## 0.00000000 0.00000000 0.00000000 0.08333333 0.07692308 0.25000000 0.08333333
## 0.9 1 1.1 1.2
## 0.12500000 0.00000000 0.00000000 0.50000000tapply(sim$IndicGen_aov, sim$SA_NonGen, mean)## 0.2 0.3 0.4 0.5 0.6 0.7 0.8
## 0.00000000 0.00000000 0.00000000 0.08333333 0.07692308 0.25000000 0.08333333
## 0.9 1 1.1 1.2
## 0.12500000 0.00000000 0.00000000 0.50000000## Local adaptation (non-genetic)
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)## 0.2 0.3 0.4 0.5 0.6 0.7 0.8
## 0.00000000 0.33333333 0.08333333 0.08333333 0.42307692 0.37500000 0.58333333
## 0.9 1 1.1 1.2
## 0.50000000 0.25000000 0.66666667 1.00000000tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)## 0.2 0.3 0.4 0.5 0.6 0.7 0.8
## 0.00000000 0.00000000 0.08333333 0.08333333 0.42307692 0.37500000 0.58333333
## 0.9 1 1.1 1.2
## 0.50000000 0.25000000 0.66666667 1.000000003.2 Poisson distribution
3.2.1 Balanced
rm(sim)
rm(data_proportion_gen, data_proportion_nongen,
sim_select_genetic, sim_select_nongenetic,
data_prop, data_prop_gen, data_prop_nongen)
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho = 0, rho_ng = 0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho = -0.05, rho_ng = -0.05,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho = -0.5, rho_ng = -0.5,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
#Add indic SA
sim <- add_indic_sign(sim)
#Add number of simul for each SA value to keep only 100 simuls per SA value
sim <- add_sim_number_SA(sim)
## Local adaptation (genetic)
tapply(sim$IndicGen[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], length)## -2.2 -2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6
## 1 1 3 3 6 8 8 15 26 38 46 79 93 100 100 100
## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
## 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6
## 100 100 100 100 100 100 100 100 100 100 75 56 41 13 11 11
## 2.7 2.9 3
## 6 1 1tapply(sim$IndicGen_aov[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], mean)## -2.2 -2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2
## 0.0000000 0.0000000 0.0100000 0.0000000 0.0000000 0.0000000 0.0000000 0.0100000
## 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## 0.0400000 0.0100000 0.0400000 0.1000000 0.0600000 0.1000000 0.1600000 0.3900000
## 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
## 0.4800000 0.5500000 0.7300000 0.7700000 0.8400000 0.9000000 0.9000000 0.9000000
## 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6
## 0.9200000 0.9400000 0.9733333 0.9821429 1.0000000 1.0000000 1.0000000 1.0000000
## 2.7 2.9 3
## 1.0000000 1.0000000 1.0000000## Local adaptation (non-genetic)
tapply(sim$IndicNonGen[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)## -2.4 -2.1 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## 0.0000000 0.0200000 0.0300000 0.0400000 0.0800000 0.0900000 0.1500000 0.3300000
## 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
## 0.3700000 0.4600000 0.5700000 0.6300000 0.6300000 0.7300000 0.7800000 0.8400000
## 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6
## 0.9100000 0.9600000 0.9387755 0.9464286 0.9722222 1.0000000 0.8571429 1.0000000
## 2.7 2.8 2.9 3
## 1.0000000 1.0000000 1.0000000 1.0000000tapply(sim$IndicNonGen_aov[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)## -2.4 -2.1 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0100000 0.0100000 0.0100000
## 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## 0.0000000 0.0400000 0.0200000 0.0400000 0.0800000 0.1100000 0.1600000 0.3100000
## 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
## 0.3900000 0.4400000 0.5400000 0.6500000 0.6400000 0.7300000 0.7900000 0.8400000
## 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6
## 0.9200000 0.9500000 0.9387755 0.9464286 0.9722222 0.9629630 0.8571429 1.0000000
## 2.7 2.8 2.9 3
## 1.0000000 1.0000000 1.0000000 1.0000000tapply(sim$IndicGen, sim$SA_Gen, length)## -2.2 -2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6
## 1 1 3 3 6 8 8 15 26 38 46 79 93 117 159 228
## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## 280 283 316 289 353 303 361 313 326 405 465 645 791 917 1058 1056
## 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6
## 1079 1062 917 723 609 474 330 287 175 137 75 56 41 13 11 11
## 2.7 2.9 3
## 6 1 1tapply(sim$IndicNonGen, sim$SA_NonGen, length)## -2.4 -2.1 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6
## 1 1 1 1 5 9 9 19 24 33 50 89 100 131 168 201
## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## 253 277 303 299 349 348 349 320 303 374 454 619 833 985 1080 1088
## 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6
## 1094 976 869 760 628 444 363 249 162 131 98 56 36 27 14 7
## 2.7 2.8 2.9 3
## 4 3 2 13.2.1.0.1 TEST
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)## -2.2 -2 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1 -0.9 -0.8 -0.7 -0.6
## 1 1 3 3 6 8 8 15 26 38 46 79 93 117 159 228
## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## 280 283 316 289 353 303 361 313 326 405 465 645 791 917 1058 1056
## 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6
## 1079 1062 917 723 609 474 330 287 175 137 75 56 41 13 11 11
## 2.7 2.9 3
## 6 1 1sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
data_prop$design <- 0
data_prop_gen <- data_prop
#Plot
plot_genetic_poisson<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Poisson data") +
theme_LO_sober +
ggtitle("Poisson data") +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_genetic_poisson
3.2.1.0.2 TEST
#################################
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 0
data_prop_nongen <- data_prop
#Plot
plot_nongenetic_poisson<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
theme_LO_sober +
ggtitle("Poisson data") +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_nongenetic_poisson
3.2.2 Unbalanced (low)
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
design = "unbalanced", disp_design = 1,
rho = 0, rho_ng = 0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
design = "unbalanced", disp_design = 1,
rho = -0.5, rho_ng = -0.5,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
design = "unbalanced", disp_design = 1,
rho = -0.05, rho_ng = -0.05,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA
sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
data_prop$design <- 1
data_prop_gen <- rbind(data_prop_gen,data_prop)
#Plot
plot_genetic_poisson_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Poisson data") +
theme_LO_sober +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_genetic_poisson_unbalanced
#################################
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 1
data_prop_nongen <- rbind(data_prop_nongen,data_prop)
#Plot
plot_nongenetic_poisson_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
theme_LO_sober +
ggtitle("Poisson data") +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_nongenetic_poisson_unbalanced
3.2.3 Unbalanced (high)
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
design = "unbalanced", disp_design = 4,
rho = 0, rho_ng = 0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
design = "unbalanced", disp_design = 4,
rho = -0.5, rho_ng = -0.5,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson",
design = "unbalanced", disp_design = 4,
rho = -0.05, rho_ng = -0.05,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA
sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
data_prop$design <- 4
data_prop_gen <- rbind(data_prop_gen,data_prop)
#Plot
plot_genetic_poisson_unbalanced_high<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Poisson data") +
theme_LO_sober +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_genetic_poisson_unbalanced_high
#################################
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 4
data_prop_nongen <- rbind(data_prop_nongen,data_prop)
#Plot
plot_nongenetic_poisson_unbalanced_high<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
theme_LO_sober +
ggtitle("Poisson data") +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_nongenetic_poisson_unbalanced_high
3.2.4 Plot
data_proportion_gen<-tidyr::gather(data_prop_gen, "method", "proportion", 2:3)
data_proportion_gen$method <- plyr::revalue(data_proportion_gen$method,
c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_gen$design <- as.factor(data_proportion_gen$design)
#Plot genetic
plot_genetic_poisson_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_gen,
aes(x = val_SA_genet, y = proportion, color = method, linetype = design)) +
geom_point(size = 2) +
geom_line(size = 1) +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Poisson data") +
theme_LO_sober +
labs(color = "Methods", linetype = "Unbalanced\ndesign") +
scale_color_manual(values = c("#459B76","#DCA237")) +
scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_genetic_poisson_balancedvsunbalanced
#Data non-genetic
data_proportion_nongen<-tidyr::gather(data_prop_nongen, "method", "proportion", 2:3)
data_proportion_nongen$method <- plyr::revalue(data_proportion_nongen$method,
c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_nongen$design <- as.factor(data_proportion_nongen$design)
#Plot non-genetic
plot_nongenetic_poisson_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_nongen,
aes(x = val_SA_nongenet, y = proportion, color = method, linetype = design)) +
geom_point(size = 2) +
geom_line(size = 1) +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
ggtitle("Poisson data") +
theme_LO_sober +
ggtitle("Poisson data") +
labs(color = "Methods", linetype = "Unbalanced\ndesign") +
scale_color_manual(values = c("#459B76","#DCA237")) +
scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_nongenetic_poisson_balancedvsunbalanced
3.2.5 False positive rate analyses
3.2.5.1 No genetic and no non-genetic effects
sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho=0, rho_ng=0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
## Local adaptation (genetic)
mean(sim$IndicGen)## [1] 0.03mean(sim$IndicGen_aov)## [1] 0.03## Local adaptation (non-genetic)
mean(sim$IndicNonGen)## [1] 0.01mean(sim$IndicNonGen_aov)## [1] 0.013.2.5.2 Genetic effects and no non-genetic effects
sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho=-0.1, rho_ng=0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
## Local adaptation (genetic)
tapply(sim$IndicGen, sim$SA_Gen, mean)## 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## 0.0000000 0.0000000 0.0000000 0.1250000 0.0000000 0.3333333 0.2500000 0.3076923
## 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
## 0.3636364 0.4444444 0.3750000 0.6666667 0.3750000 0.6000000 1.0000000 1.0000000
## 1.9 2
## 1.0000000 1.0000000tapply(sim$IndicGen_aov, sim$SA_Gen, mean)## 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## 0.0000000 0.0000000 0.0000000 0.1250000 0.0000000 0.3333333 0.2500000 0.3076923
## 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
## 0.3636364 0.4444444 0.3750000 0.6666667 0.3750000 0.6000000 1.0000000 1.0000000
## 1.9 2
## 1.0000000 1.0000000## Local adaptation (non-genetic)
mean(sim$IndicNonGen)## [1] 0.04mean(sim$IndicNonGen_aov)## [1] 0.043.2.5.3 Non-genetic effects and no genetic effects
sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "poisson",design = "balanced",
rho=0, rho_ng=-0.1,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
## Local adaptation (genetic)
mean(sim$IndicGen)## [1] 0.12mean(sim$IndicGen_aov)## [1] 0.12## Local adaptation (non-genetic)
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)## 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## 0.0000000 0.0000000 0.0000000 0.2500000 0.1250000 0.0000000 0.4666667 0.4615385
## 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
## 0.2000000 0.5000000 0.5000000 0.1666667 0.6666667 0.0000000 0.6666667 1.0000000
## 1.9 2
## 1.0000000 1.0000000tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)## 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
## 0.0000000 0.0000000 0.0000000 0.2500000 0.1250000 0.0000000 0.4000000 0.4615385
## 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8
## 0.2000000 0.4000000 0.5000000 0.1666667 0.6666667 0.0000000 1.0000000 1.0000000
## 1.9 2
## 1.0000000 1.00000003.3 Binomial distribution
3.3.1 Balanced
rm(sim)
rm(data_proportion_gen, data_proportion_nongen,
sim_select_genetic, sim_select_nongenetic,
data_prop, data_prop_gen, data_prop_nongen)
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial", design = "balanced",
rho = -0.05, rho_ng = -0.05,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_databinomial,
sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial,
sigma = sigma_databinomial)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial", design = "balanced",
rho = 0, rho_ng = 0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_databinomial,
sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial,
sigma = sigma_databinomial)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial", design = "balanced",
rho = -0.5, rho_ng = -0.5,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_databinomial,
sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial,
sigma = sigma_databinomial)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
#Add indic SA
sim <- add_indic_sign(sim)
#Add number of simul for each SA value to keep only 100 simuls per SA value
sim <- add_sim_number_SA(sim)
## Local adaptation (genetic)
tapply(sim$IndicGen[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], length)## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6
## 1 11 79 100 100 100 100 100 100 100 100 7tapply(sim$IndicGen_aov[sim$index_SA_Gen<101], sim$SA_Gen[sim$index_SA_Gen<101], mean)## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0300000 0.4500000 0.6100000
## 0.3 0.4 0.5 0.6
## 0.6900000 0.7600000 0.8100000 0.8571429## Local adaptation (non-genetic)
tapply(sim$IndicNonGen[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6
## 0.00 0.00 0.00 0.00 0.00 0.00 0.04 0.24 0.33 0.61 0.74 0.75tapply(sim$IndicNonGen_aov[sim$index_SA_NonGen<101], sim$SA_NonGen[sim$index_SA_NonGen<101], mean)## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2
## 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0100000 0.0500000 0.3000000
## 0.3 0.4 0.5 0.6
## 0.3300000 0.6200000 0.6500000 0.8333333tapply(sim$IndicGen, sim$SA_Gen, length)## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6
## 1 11 79 409 1245 1550 2351 4904 3376 943 124 7tapply(sim$IndicNonGen, sim$SA_NonGen, length)## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6
## 1 6 75 449 1173 1601 2281 5070 3254 939 139 12####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)## -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6
## 1 11 79 409 1245 1550 2351 4904 3376 943 124 7sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
data_prop$design <- 0
data_prop_gen <- data_prop
#Plot
plot_genetic_binomial<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Binomial data") +
theme_LO_sober +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_genetic_binomial
#################################
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 0
data_prop_nongen <- data_prop
#Plot
plot_nongenetic_binomial<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Balanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Balanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
theme_LO_sober +
ggtitle("Binomial data") +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_nongenetic_binomial
3.3.2 Unbalanced (low)
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
design = "unbalanced", disp_design = 1,
rho = -0.5, rho_ng = -0.5,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_databinomial,
sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial,
sigma = sigma_databinomial)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
design = "unbalanced", disp_design = 1,
rho = 0, rho_ng = 0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_databinomial,
sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial,
sigma = sigma_databinomial)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
design = "unbalanced", disp_design = 1,
rho = -0.05, rho_ng = -0.05,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_databinomial,
sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial,
sigma = sigma_databinomial)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA
sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
data_prop$design <- 1
data_prop_gen <- rbind(data_prop_gen,data_prop)
#Plot
plot_genetic_binomial_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Binomial data") +
theme_LO_sober +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_genetic_binomial_unbalanced
#################################
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 1
data_prop_nongen <- rbind(data_prop_nongen,data_prop)
#Plot
plot_nongenetic_binomial_unbalanced<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
theme_LO_sober +
ggtitle("Binomial data") +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_nongenetic_binomial_unbalanced
3.3.3 Unbalanced (high)
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
design = "unbalanced", disp_design = 4,
rho = 0, rho_ng = 0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_databinomial,
sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial,
sigma = sigma_databinomial)
sim1 <- as.data.frame(t(sim1))
#Perform siluations
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
design = "unbalanced", disp_design = 4,
rho = -0.05, rho_ng = -0.05,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_databinomial,
sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial,
sigma = sigma_databinomial)
sim2 <- as.data.frame(t(sim2))
#Perform siluations
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "binomial",
design = "unbalanced", disp_design = 4,
rho = -0.5, rho_ng = -0.5,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_databinomial,
sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial,
sigma = sigma_databinomial)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
# tapply(sim$index_SA_Gen, sim$SA_Gen, length)
# nb_per_SA
sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
data_prop$design <- 4
data_prop_gen <- rbind(data_prop_gen,data_prop)
#Plot
plot_genetic_binomial_unbalanced_high<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_genet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Binomial data") +
theme_LO_sober +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_genetic_binomial_unbalanced_high
#################################
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
data_prop$design <- 4
data_prop_nongen <- rbind(data_prop_nongen,data_prop)
#Plot
plot_nongenetic_binomial_unbalanced_high<-ggplot2::ggplot(data = data_prop, aes(x=val_SA_nongenet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 4) +
geom_line(aes(y=prop_val_sign, color="LMER", linetype = "Unbalanced"), size = 1.5) +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova", linetype = "Unbalanced"), size = 0.75) +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
theme_LO_sober +
ggtitle("Binomial data") +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_nongenetic_binomial_unbalanced_high
3.3.4 Plot
data_proportion_gen<-tidyr::gather(data_prop_gen, "method", "proportion", 2:3)
data_proportion_gen$method <- plyr::revalue(data_proportion_gen$method,
c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_gen$design <- as.factor(data_proportion_gen$design)
#Plot genetic
plot_genetic_binomial_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_gen,
aes(x = val_SA_genet, y = proportion, color = method, linetype = design)) +
geom_point(size = 2) +
geom_line(size = 1) +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Binomial data") +
theme_LO_sober +
labs(color = "Methods", linetype = "Unbalanced\ndesign") +
scale_color_manual(values = c("#459B76","#DCA237")) +
scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_genetic_binomial_balancedvsunbalanced
#Data non-genetic
data_proportion_nongen<-tidyr::gather(data_prop_nongen, "method", "proportion", 2:3)
data_proportion_nongen$method <- plyr::revalue(data_proportion_nongen$method,
c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_nongen$design <- as.factor(data_proportion_nongen$design)
#Plot non-genetic
plot_nongenetic_binomial_balancedvsunbalanced<-ggplot2::ggplot(data = data_proportion_nongen,
aes(x = val_SA_nongenet, y = proportion, color = method, linetype = design)) +
geom_point(size = 2) +
geom_line(size = 1) +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
ggtitle("Binomial data") +
theme_LO_sober +
labs(color = "Methods", linetype = "Unbalanced\ndesign") +
scale_color_manual(values = c("#459B76","#DCA237")) +
scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_nongenetic_binomial_balancedvsunbalanced
3.3.5 False positive rate analyses
3.3.5.1 No genetic and no non-genetic effects
sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "binomial", design = "balanced",
rho=0, rho_ng=0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_databinomial,
sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial,
sigma = sigma_databinomial)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
## Local adaptation (genetic)
mean(sim$IndicGen)## [1] 0.02mean(sim$IndicGen_aov)## [1] 0.02## Local adaptation (non-genetic)
mean(sim$IndicNonGen)## [1] 0mean(sim$IndicNonGen_aov)## [1] 0.023.3.5.2 Genetic effects and no non-genetic effects
sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "binomial", design = "balanced",
rho=-0.1, rho_ng=0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_databinomial,
sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial,
sigma = sigma_databinomial)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
## Local adaptation (genetic)
tapply(sim$IndicGen, sim$SA_Gen, mean)## 0.1 0.2 0.3 0.4
## 0.0500000 0.2926829 0.4062500 0.5714286tapply(sim$IndicGen_aov, sim$SA_Gen, mean)## 0.1 0.2 0.3 0.4
## 0.0500000 0.2926829 0.4062500 0.5714286## Local adaptation (non-genetic)
mean(sim$IndicNonGen)## [1] 0mean(sim$IndicNonGen_aov)## [1] 03.3.5.3 non-genetic effects and no genetic effects
###
sim <- sapply(1:nb_simul_fpr, simul_fitnessdata, distrib = "binomial", design = "balanced",
rho=0, rho_ng=-0.1,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_databinomial,
sdfruithab = sdfruithab_databinomial, sdfruithab_ng = sdfruithab_ng_databinomial,
sigma = sigma_databinomial)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
## Local adaptation (genetic)
mean(sim$IndicGen)## [1] 0.07mean(sim$IndicGen_aov)## [1] 0.07## Local adaptation (non-genetic)
tapply(sim$IndicNonGen, sim$SA_NonGen, mean)## 0.1 0.2 0.3 0.4
## 0.0000000 0.2075472 0.3333333 0.6000000tapply(sim$IndicNonGen_aov, sim$SA_NonGen, mean)## 0.1 0.2 0.3 0.4
## 0.0000000 0.3018868 0.4074074 0.60000004 PLOT Balanced vs Unbalanced
POWER_ALL<-cowplot::plot_grid(plot_genetic_normal_balancedvsunbalanced,
plot_genetic_poisson_balancedvsunbalanced,
plot_genetic_binomial_balancedvsunbalanced,
plot_nongenetic_normal_balancedvsunbalanced,
plot_nongenetic_poisson_balancedvsunbalanced,
plot_nongenetic_binomial_balancedvsunbalanced,
ncol = 3, nrow = 2,
scale = c(1, 1))
POWER_ALL
#Save plot
name_plot<-paste0("Simul_powertest_ALLDATA_BALANCED_and_UNBALANCED",
nb_simul,"nbsimul_",nb_per_SA,
"nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
# POWER_ALL, base_height = 20/cm(1),
# base_width = 40/cm(1), dpi = 1200)5 BOX EFFECT Absent vs Present
5.1 Non box
rm(sim)
rm(data_proportion_gen, data_proportion_nongen,
sim_select_genetic, sim_select_nongenetic,
data_prop, data_prop_gen, data_prop_nongen)
#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho = -0.5, rho_ng = -0.5,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho = -0.05, rho_ng = -0.05,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho = 0, rho_ng = 0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)
colors <- c("LMER" = "#DCA237", "Anova" = "#459B76")
lines_box <- c("Present" = "dotted", "Absent" = "solid")
####################################
###### Plot genetic
####################################
sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
data_prop$sd_box <- 0
data_prop_gen <- data_prop
#################################
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
data_prop$sd_box <- 0
data_prop_nongen <- data_prop5.2 Box (low)
#With lower value of sdbox (close to sdpop):
sdbox_datapoisson_box## BoxID
## 0.4464253sigma_datapoisson## [1] 0.8727978sdpop_datapoisson## Population
## 0.504277sdfruithab_datapoisson## Original_environment:Test_environment:IndicG0
## 0.8111606#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho = -0.5, rho_ng = -0.5,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 0.5,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho = 0, rho_ng = 0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 0.5,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho = -0.05, rho_ng = -0.05,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = sdbox_datapoisson_box,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)
####################################
###### Plot genetic
####################################
sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
data_prop$sd_box <- 0.5
data_prop_gen <- rbind(data_prop_gen,data_prop)
#################################
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
data_prop$sd_box <- 0.5
data_prop_nongen <- rbind(data_prop_nongen,data_prop)5.3 Box (high)
#With lower value of sdbox (close to sdpop):
sdbox_datapoisson_box## BoxID
## 0.4464253sigma_datapoisson## [1] 0.8727978sdpop_datapoisson## Population
## 0.504277sdfruithab_datapoisson## Original_environment:Test_environment:IndicG0
## 0.8111606#Perform siluations
sim1 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho = -0.5, rho_ng = -0.5,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 1.5,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim1 <- as.data.frame(t(sim1))
sim2 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho = 0, rho_ng = 0,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 1.5,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim2 <- as.data.frame(t(sim2))
sim3 <- sapply(1:nb_simul, simul_fitnessdata, distrib = "poisson", design = "balanced",
rho = -0.05, rho_ng = -0.05,
npop_per_fruit = npop_per_fruit_data, nfruit = nfruit_data, nhab = nhab_data,
nrep = nrep_data, sdpop = sdpop_datapoisson, sdbox = 1.5,
sdfruithab = sdfruithab_datapoisson, sdfruithab_ng = sdfruithab_ng_datapoisson,
sigma = sigma_datapoisson)
sim3 <- as.data.frame(t(sim3))
# Concatenate value
sim <- rbind(sim1,sim2, sim3)
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)
####################################
###### Plot genetic
####################################
sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
data_prop$sd_box <- 1.5
data_prop_gen <- rbind(data_prop_gen,data_prop)
#################################
### Plot non-genetic
#################################
#Save the same number of simuls for each SA value
sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
data_prop$sd_box <- 1.5
data_prop_nongen <- rbind(data_prop_nongen,data_prop)5.4 Plot
#################################
### Common plot
#################################
data_proportion_gen<-tidyr::gather(data_prop_gen, "method", "proportion", 2:3)
data_proportion_gen$method <- plyr::revalue(data_proportion_gen$method,
c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_gen$sd_box <- as.factor(data_proportion_gen$sd_box)
#Plot genetic
plot_genetic_poisson_ALL_box<-ggplot2::ggplot(data = data_proportion_gen,
aes(x = val_SA_genet, y = proportion, color = method, linetype = sd_box)) +
geom_point(size = 2) +
geom_line(size = 1, alpha = 0.7) +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Poisson data") +
theme_LO_sober +
labs(color = "Methods", linetype = "Box effect\n(sd box)") +
scale_color_manual(values = c("#459B76","#DCA237")) +
scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_genetic_poisson_ALL_box
#Data non-genetic
data_proportion_nongen<-tidyr::gather(data_prop_nongen, "method", "proportion", 2:3)
data_proportion_nongen$method <- plyr::revalue(data_proportion_nongen$method,
c("prop_val_sign"="LMER", "prop_val_sign_aov"="Anova"))
data_proportion_nongen$sd_box <- as.factor(data_proportion_nongen$sd_box)
#Plot non-genetic
plot_nongenetic_poisson_ALL_box<-ggplot2::ggplot(data = data_proportion_nongen,
aes(x = val_SA_nongenet, y = proportion, color = method, linetype = sd_box)) +
geom_point(size = 2) +
geom_line(size = 1) +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
ggtitle("Poisson data") +
theme_LO_sober +
labs(color = "Methods", linetype = "Box effect\n(sd box)") +
scale_color_manual(values = c("#459B76","#DCA237")) +
scale_linetype_manual(values = c("solid","dashed","dotted"))
plot_nongenetic_poisson_ALL_box
power_poisson_box<-cowplot::plot_grid(plot_genetic_poisson_ALL_box,
plot_nongenetic_poisson_ALL_box,
labels=c("A", "B"),
label_size = 15,
ncol =1, nrow = 2,
# hjust = 0, vjust = 1,
scale = c(1, 1))
power_poisson_box
#Save plot
name_plot<-paste0("Simul_powertest_poisson_boxeffect_",
nb_simul,"nbsimul_",nb_per_SA,
"nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
# power_poisson_box, base_height = 20/cm(1), base_width = 14/cm(1), dpi = 1200)6 SIMUL WITH REAL UNBALANCED DESIGN
6.1 Dataset
#unbalanced dataset
data_PERF <- read.table(file=here::here("data", "DATACOMPLET_PERF.csv"), sep=",", header=TRUE)
data_PERF$Original_environment<-as.factor(data_PERF$Original_environment)
data_PERF <- data_PERF[data_PERF$Original_environment!="WT3",]
data_PERF <- data_PERF[data_PERF$Test_environment!="GF",]
data_PERF <- data_PERF[data_PERF$Test_environment!="Grape",] #mvrnom doesn't work with unbalanced matrix
data_PERF <- droplevels(data_PERF)
tapply(data_PERF$Population, list(data_PERF$Generation, data_PERF$Population), length)## Blackberry31 Blackberry32 Blackberry33 Blackberry34 Blackberry35
## G0 38 42 25 4 NA
## G2 27 39 44 11 46
## Blackberry36 Blackberry37 Blackberry38 Blackberry39 Blackberry40
## G0 8 15 21 10 15
## G2 17 33 5 24 30
## Blackberry43 Blackberry44 Blackberry45 Cherry103 Cherry104 Cherry3 Cherry47
## G0 8 12 43 3 39 9 36
## G2 10 30 4 19 16 9 39
## Cherry50 Cherry51 Cherry52 Cherry6 Cherry7 Strawberry42 Strawberry44
## G0 36 48 45 36 6 15 21
## G2 11 1 4 7 37 46 30
## Strawberry53
## G0 6
## G2 636.2 Normal distribution
####################################
###### Simulation
####################################
#Simul data
sim <- sapply(1:nb_simul, simul_fitnessdata_unbalanced, distrib = "normal",
unbalanced_dataset = data_PERF,
sdpop = 0.5, sdfruithab = 0.5, sdfruithab_ng = 0.5,
rho = -0.05, rho_ng = -0.05,
sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)## 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4
## 2 58 240 568 819 944 827 675 421 240 127 50 21 8#Save the same number of simuls for each SA value
nb_per_SA = 50
sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
#Plot
plot_genetic_normal_unbalanced_real<-ggplot2::ggplot(data = data_prop,
aes(x=val_SA_genet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) +
geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Normal data (unbalanced experimental design)") +
theme_LO_sober +
labs(color = "Methods", linetype = "Design") +
scale_color_manual(values = colors) +
scale_linetype_manual(values = lines)
plot_genetic_normal_unbalanced_real
#################################
### Plot non-genetic
#################################
#Number of simuls for each SA non-genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)## 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5
## 1 60 236 528 890 969 826 617 430 215 143 47 26 8 4#Save the same number of simuls for each SA value
nb_per_SA ## [1] 50sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
#Plot
plot_nongenetic_normal_unbalanced_real<-ggplot2::ggplot(data = data_prop,
aes(x=val_SA_nongenet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) +
geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
theme_LO_sober +
ggtitle("Normal data (unbalanced experimental design)") +
labs(color = "Methods") +
scale_color_manual(values = colors)
plot_nongenetic_normal_unbalanced_real
6.3 Poisson distribution
####################################
###### Simulation
####################################
#Simul data
#nb_simul <- 1000
sim <- sapply(1:nb_simul, simul_fitnessdata_unbalanced, distrib = "poisson",
unbalanced_dataset = data_PERF, sdpop = 0.5,
sdfruithab = 1, sdfruithab_ng = 1, rho = -0.05, rho_ng = -0.05, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)## 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2
## 3 26 66 123 195 289 353 422 461 466 451 394 387 347 253 233 166 109 86 71
## 2.3 2.4 2.5 2.6 2.7 2.8 2.9
## 34 29 15 9 6 4 2#Save the same number of simuls for each SA value
# nb_per_SA=10
sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
#Plot
plot_genetic_poisson_unbalanced_real<-ggplot2::ggplot(data = data_prop,
aes(x=val_SA_genet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) +
geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Poisson data (unbalanced experimental design)") +
theme_LO_sober +
labs(color = "Methods") +
scale_color_manual(values = colors)
plot_genetic_poisson_unbalanced_real
#################################
### Plot non-genetic
#################################
#Number of simuls for each SA non-genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)## 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1
## 1 4 29 61 123 174 256 386 449 495 478 443 405 379 326 239 201 174 107 93
## 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3
## 63 40 28 17 14 7 2 3 3#Save the same number of simuls for each SA value
nb_per_SA ## [1] 50sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
colors<-c("LMER" = "#DCA237", "Anova" = "#459B76")
#Plot
plot_nongenetic_poisson_unbalanced_real<-ggplot2::ggplot(data = data_prop,
aes(x=val_SA_nongenet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) +
geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
theme_LO_sober +
ggtitle("Poisson data (unbalanced experimental design)") +
labs(color = "Methods") +
scale_color_manual(values = colors)
plot_nongenetic_poisson_unbalanced_real
6.4 Binomial distribution
####################################
###### Simulation
####################################
#Simul data
#nb_simul <- 1000
sim <- sapply(1:nb_simul, simul_fitnessdata_unbalanced, distrib = "binomial",
unbalanced_dataset = data_PERF, sdpop = 0.5,
sdfruithab = 0.5, sdfruithab_ng = 0.5,
rho = -0.05, rho_ng = -0.05, sigma = 0.5)
sim <- as.data.frame(t(sim))
sim <- add_indic_sign(sim)
sim <- add_sim_number_SA(sim)
####################################
###### Plot genetic
####################################
#Number of simuls for each SA value
tapply(sim$index_SA_Gen, sim$SA_Gen, length)## 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4
## 2 58 240 568 819 944 827 675 421 240 127 50 21 8#Save the same number of simuls for each SA value
# nb_per_SA=10
sim_select_genetic <- select_sample_simul(dataset_sim = sim,
type_SA = "genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_genet = sort(unique(sim_select_genetic$SA_Gen)),
prop_val_sign = tapply(sim_select_genetic$IndicGen,
sim_select_genetic$SA_Gen, mean),
prop_val_sign_aov = tapply(sim_select_genetic$IndicGen_aov,
sim_select_genetic$SA_Gen, mean))
#Plot
plot_genetic_Binomial_unbalanced_real<-ggplot2::ggplot(data = data_prop,
aes(x=val_SA_genet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) +
geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") +
ylab("Proportion of simulations with\nsignificant genetic local adaption") +
xlab("Intensity of genetic local adaptation (SA)") +
ggtitle("Binomial data (unbalanced experimental design)") +
theme_LO_sober +
labs(color = "Methods") +
scale_color_manual(values = colors)
plot_genetic_Binomial_unbalanced_real
#################################
### Plot non-genetic
#################################
#Number of simuls for each SA non-genetic value
tapply(sim$index_SA_NonGen, sim$SA_NonGen, length)## 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5
## 1 60 236 528 890 969 826 617 430 215 143 47 26 8 4#Save the same number of simuls for each SA value
nb_per_SA ## [1] 50sim_select_nongenetic <- select_sample_simul(dataset_sim = sim,
type_SA = "non-genetic",
nb_simul_per_SA = nb_per_SA,
minimum_SA = "0")
#Compute proportion of significant simul for each SA values
data_prop <- data.frame(val_SA_nongenet = sort(unique(sim_select_nongenetic$SA_NonGen)),
prop_val_sign = tapply(sim_select_nongenetic$IndicNonGen,
sim_select_nongenetic$SA_NonGen, mean),
prop_val_sign_aov = tapply(sim_select_nongenetic$IndicNonGen_aov,
sim_select_nongenetic$SA_NonGen, mean))
colors<-c("LMER" = "#DCA237", "Anova" = "#459B76")
#Plot
plot_nongenetic_Binomial_unbalanced_real<-ggplot2::ggplot(data = data_prop,
aes(x=val_SA_nongenet)) +
geom_point(aes(y=prop_val_sign, color="LMER"), size = 2.5) +
geom_line(aes(y=prop_val_sign, color="LMER"), size = 1, linetype = "dotted") +
geom_point(aes(y=prop_val_sign_aov, color="Anova"), size = 2.5) +
geom_line(aes(y=prop_val_sign_aov, color="Anova"), size = 1, linetype = "dotted") +
ylab("Proportion of simulations with\nsignificant non-genetic local adaption") +
xlab("Intensity of non-genetic local adaptation (SA)") +
theme_LO_sober +
ggtitle("Binomial data (unbalanced experimental design)") +
labs(color = "Methods") +
scale_color_manual(values = colors)
plot_nongenetic_Binomial_unbalanced_real
6.5 All plots unbalanced
power_plot_all_unbalanced_real<-cowplot::plot_grid(plot_nongenetic_normal_unbalanced_real,
plot_nongenetic_poisson_unbalanced_real,
plot_nongenetic_Binomial_unbalanced_real,
plot_genetic_normal_unbalanced_real,
plot_genetic_poisson_unbalanced_real,
plot_genetic_Binomial_unbalanced_real,
#labels=c("A", "B"),
#label_size = 15,
ncol =3, nrow = 2,
# hjust = 0, vjust = 1,
scale = c(1, 1))
power_plot_all_unbalanced_real
#Save plot
name_plot<-paste0("Simul_powertest_ALLDATA_UNBALANCED_real",
nb_simul,"nbsimul_",nb_per_SA,
"nbsimulperSA.pdf")
# cowplot::save_plot(file =here::here("figures",name_plot ),
# power_plot_all_unbalanced_real, base_height = 20/cm(1),
# base_width = 42/cm(1), dpi = 1200)
#